Nesta lista, vamos trabalhar com o seguinte processo estocástico AR(1):
\[ z_t = \rho z_{t-1} + \varepsilon_t , \]
com \(\varepsilon_t \sim N(\mu,\sigma^2)\). Inicialmente, assumiremos que \(\mu=0\), \(\rho = 0.95\) e \(\sigma = 0.007\), calibração de Cooley & Prescott (1995).
Ao discretizar um processo que assume valores contínuos, é necessário determinar o intervalo do grid, o número de pontos, e sua localização. Além disso, é necessário determinar probabiliades condicionais de transição entre os estados de forma consistente com o processo original.
Pelo método de Tauchen, o intervalo de um grid de \(N\) pontos é determinado por \([\mu - \theta_N, \mu + \theta_N]\), em que \(\theta_N = m \frac{\sigma}{\sqrt{1-\rho^2}}\) é o desvio padrão do processo estocástico, multiplicado por um parâmetro de escala \(m\). Os demais \(N-2\) pontos são distribuídos uniformemente no intervalo.
# Parâmetros Iniciais
N = 9
rho = 0.95
sigma = 0.007
mu = 0
m = 3
tauchen_grid = function(N, m, sigma, rho) {
thetaN = m * sigma / sqrt(1-rho^2)
seq(-thetaN, thetaN, length.out = N)
}
thetaT = tauchen_grid(N, m, sigma, rho)
as.matrix(thetaT)## [,1]
## [1,] -0.06725382
## [2,] -0.05044037
## [3,] -0.03362691
## [4,] -0.01681346
## [5,] 0.00000000
## [6,] 0.01681346
## [7,] 0.03362691
## [8,] 0.05044037
## [9,] 0.06725382Seja \(\theta_{i,t}\) um estado no período \(t\), e \(\theta_{j, t+1}\) um estado no período \(t+1\). No método de Tauchen, as probabilidades de transição dos pontos interiores do grid são dadas por
\[ p_{i,j} = \Phi \left( \frac{\theta_j + \Delta \theta/2 - (1-\rho) \mu - \rho \theta_i}{\sigma} \right) - \Phi \left(\frac{\theta_j - \Delta \theta/2 - (1-\rho) \mu - \rho \theta_i}{\sigma}\right)\]
em que \(\Phi\) é a CDF da \(N(\mu, \sigma^2)\). Nos cantos do grid, temos que
\[ p_{i,1} = \Phi \left(\frac{\theta_1 - (1-\rho) \mu - \rho \theta_i + \Delta \theta/2}{\sigma} \right) \]
\[ p_{i,N} = 1 - \Phi \left(\frac{\theta_N - (1-\rho) \mu - \rho \theta_i - \Delta \theta/2}{\sigma} \right).\]
tauchen_P = function(grid, rho, sigma, mu) {
N = length(grid)
delta = (max(grid) - min(grid)) / (N-1)
PT = matrix(NA_real_, N, N)
for(i in 1:N) {
for(j in 1:N) {
if(j == 1) {
PT[i,j] = pnorm( (grid[1] - (1-rho)*mu - rho*grid[i] + delta/2) / sigma )
} else if(j==N) {
PT[i,j] = 1 - pnorm( (grid[N] - (1-rho)*mu - rho*grid[i] - delta/2) / sigma )
} else {
PT[i,j] =
pnorm((grid[j] + delta/2 - (1-rho)*mu - rho*grid[i]) / sigma) -
pnorm((grid[j] - delta/2 - (1-rho)*mu - rho*grid[i]) / sigma)
}
}
}
return(PT)
}
PT = tauchen_P(thetaT, rho, sigma, mu)
round(PT, 4)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 0.7644 0.2347 0.0009 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## [2,] 0.0592 0.7405 0.1997 0.0006 0.0000 0.0000 0.0000 0.0000 0.0000
## [3,] 0.0001 0.0747 0.7569 0.1679 0.0004 0.0000 0.0000 0.0000 0.0000
## [4,] 0.0000 0.0001 0.0931 0.7669 0.1396 0.0002 0.0000 0.0000 0.0000
## [5,] 0.0000 0.0000 0.0002 0.1147 0.7702 0.1147 0.0002 0.0000 0.0000
## [6,] 0.0000 0.0000 0.0000 0.0002 0.1396 0.7669 0.0931 0.0001 0.0000
## [7,] 0.0000 0.0000 0.0000 0.0000 0.0004 0.1679 0.7569 0.0747 0.0001
## [8,] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0006 0.1997 0.7405 0.0592
## [9,] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0009 0.2347 0.7644rowSums(PT) # Linhas da matriz devem somar 1.## [1] 1 1 1 1 1 1 1 1 1A definição do grid no método de Rouwenhorst é um caso particular do método de Tauchen, com \(m = \sqrt{N-1}\).
thetaR = tauchen_grid(N, m = sqrt(N-1), sigma, rho)
thetaR## [1] -0.06340751 -0.04755564 -0.03170376 -0.01585188 0.00000000 0.01585188
## [7] 0.03170376 0.04755564 0.06340751A matriz de transição \(P\), porém, é obtida recursivamente.
rouwen_P = function(N, rho) {
p = (1+rho)/2
P_init = matrix(c(p, 1-p, 1-p, p), 2, 2, byrow = TRUE)
for(i in 3:N) {
PR =
rbind(cbind(p*P_init, 0), 0) +
rbind(cbind(0, (1-p)*P_init), 0) +
rbind(0, cbind((1-p)*P_init, 0)) +
rbind(0, cbind(0, p*P_init))
PR = PR / rowSums(PR) # normalização
P_init = PR
}
return(PR)
}
PR = rouwen_P(N, rho)
round(PR, 4)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 0.8167 0.1675 0.0150 0.0008 0.0000 0.0000 0.0000 0.0000 0.0000
## [2,] 0.0209 0.8204 0.1469 0.0113 0.0005 0.0000 0.0000 0.0000 0.0000
## [3,] 0.0005 0.0420 0.8231 0.1261 0.0081 0.0003 0.0000 0.0000 0.0000
## [4,] 0.0000 0.0016 0.0630 0.8247 0.1051 0.0054 0.0001 0.0000 0.0000
## [5,] 0.0000 0.0001 0.0032 0.0841 0.8253 0.0841 0.0032 0.0001 0.0000
## [6,] 0.0000 0.0000 0.0001 0.0054 0.1051 0.8247 0.0630 0.0016 0.0000
## [7,] 0.0000 0.0000 0.0000 0.0003 0.0081 0.1261 0.8231 0.0420 0.0005
## [8,] 0.0000 0.0000 0.0000 0.0000 0.0005 0.0113 0.1469 0.8204 0.0209
## [9,] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0008 0.0150 0.1675 0.8167set.seed(96452)
time_span = 10000
eps = rnorm(time_span, mu, sigma)
# Processo contínuo
zc = rep(0, time_span)
for(t in 2:time_span)
zc[t] = rho * zc[t-1] + eps[t]
# Funcão para gerar processos
gen_process = function(grid, P, eps, mu, sigma) {
time_span = length(eps)
z = rep(0, time_span)
for(t in 2:time_span) {
i = which(grid == z[t-1])
j = sum(cumsum(P[i, ]) < pnorm(eps[t], mu, sigma)) + 1
z[t] = grid[j]
}
return(z)
}
zt = gen_process(thetaT, PT, eps, mu, sigma)
zr = gen_process(thetaR, PR, eps, mu, sigma)Os gráficos a seguir exibem o processo contínuo e as discretizações de Tauchen e de Rouwenhorst para 10.000 períodos e para os 500 primeiros períodos:
Visualmente, a discretização de Tauchen funciona melhor com os parâmetros dados. A discretização de Tauchen captura melhor a amplitude do processo sem perder precisão. Ambos apresentam, contudo um RMSE semelhante em relação ao processo contínuo:
RMSE = function(x,y) sqrt(mean((x-y)^2))
cat("RMSE:\n")
cbind("Tauchen" = RMSE(zc,zt), "Rouwenhorst" = RMSE(zc,zr))## RMSE:
## Tauchen Rouwenhorst
## [1,] 0.01514628 0.01515141fit_zt = lm(zt ~ shift(zt) + 0)
fit_zr = lm(zr ~ shift(zr) + 0)
coefs = c(fit_zt$coefficients,
fit_zr$coefficients)
sigmas = c(summary(fit_zt)$sigma,
summary(fit_zr)$sigma)
tbl = data.frame(method = c("Tauchen 9 pts",
"Rouwenhorst 9 pts"),
rho = coefs, sigma = sigmas)
tblPelo método de Tauchen com 9 pontos no grid, estimamos \(\hat{\sigma} = 0.008\), valor relativamente distante do \(\sigma\) verdadeiro. Aumentando o número de pontos no grid para 51, obtemos uma estimativa semelhante à obtida para o método de Rouwenhorst com 9 pontos.
thetaT = tauchen_grid(N=51, m=3, sigma, rho)
PT = tauchen_P(thetaT, rho, sigma, mu)
ztl = gen_process(thetaT, PT, eps, mu, sigma)
fit_ztl = lm(ztl ~ shift(ztl) + 0)
rbind(tbl, cbind(method = "Tauchen 25 pts", rho = fit_ztl$coefficients, sigma = summary(fit_ztl)$sigma))Vamos refazer os exercícios acima com \(\rho = 0.99\).
O método de Rouwenhorst funciona melhor para \(\rho = 0.99\) com 9 pontos no grid. Obtemos os seguintes valores de \(\hat{\rho}\) e \(\hat{\sigma}\):